/*  ----------------------------------------------------------------------------
 *  Arello Mobile
 *  Fast 3d Engine
 *  ----------------------------------------------------------------------------
 * Except where otherwise noted, this work is licensed under
 * http://creativecommons.org/licenses/by/3.0
 */

#ifndef __MINIBALL_H__
#define __MINIBALL_H__

#include <list>

namespace liba
{
namespace math
{

template <int d>
class Wrapped_array
{
private:
	double coord [d];
	
public:
	// default
	Wrapped_array()
	{}
	
	// copy from Wrapped_array
	Wrapped_array (const Wrapped_array& p)
	{
		for (int i=0; i<d; ++i)
			coord[i] = p.coord[i];
	}
	
	// copy from double*
	Wrapped_array (const double* p)
	{
		for (int i=0; i<d; ++i)
			coord[i] = p[i];
	}
	
	// assignment
	Wrapped_array& operator = (const Wrapped_array& p)
	{
		for (int i=0; i<d; ++i)
			coord[i] = p.coord[i];
		return *this;
	}
	
	// coordinate access
	double& operator [] (int i)
	{
		return coord[i];
	}
	const double& operator [] (int i) const
	{
		return coord[i];
	}
	const double* begin() const
	{
		return coord;
	}
	const double* end() const
	{
		return coord+d;
	}
};

template <int d>
class Basis
{
private:
	// types
	typedef Wrapped_array<d>            Point;

	// data members
	int                 m, s;   // size and number of support points
	double              q0[d];

	double              z[d+1];
	double              f[d+1];
	double              v[d+1][d];
	double              a[d+1][d];

	double              c[d+1][d];
	double              sqr_r[d+1];

	double*             current_c;      // points to some c[j]
	double              current_sqr_r;

	double              sqr (double r) const {return r*r;}

public:
	Basis();

	// access
	const double*       center() const;
	double              squared_radius() const;
	int                 size() const;
	int                 support_size() const;
	double              excess (const Point& p) const;

	// modification
	void                reset(); // generates empty sphere with m=s=0
	bool                push (const Point& p);
	void                pop ();

	// checking
	double              slack() const;
};

template <int d>
class Miniball
{
public:
	// types
	typedef Wrapped_array<d>                            Point;
	typedef typename std::list<Point>::iterator         It;
	typedef typename std::list<Point>::const_iterator   Cit;

private:
	// data members
	std::list<Point> L;         // STL list keeping the points
	Basis<d>    B;              // basis keeping the current ball
	It          support_end;    // past-the-end iterator of support set

	// private methods
	void        mtf_mb (It k);
	void        pivot_mb (It k);
	void        move_to_front (It j);
	double      max_excess (It t, It i, It& pivot) const;
	double      abs (double r) const {return (r>0)? r: (-r);}
	double      sqr (double r) const {return r*r;}

public:
	// construction
	Miniball() {}
	void        check_in (const Point& p);
	void        build (bool pivoting = true);

	// access
	Point       center() const;
	double      squared_radius () const;
	int         nr_points () const;
	Cit         points_begin () const;
	Cit         points_end () const;
	int         nr_support_points () const;
	Cit         support_points_begin () const;
	Cit         support_points_end () const;

	// checking
	double      accuracy (double& slack) const;
	bool        is_valid (double tolerance = 1e-15) const;
};

#include "miniball.inc"

}  // namespace math
}  // namespace liba
using namespace liba;

#endif //__MINIBALL_H__